%{
AUTHOR: Felipe Arteaga
-------------------------------------------------------------------------
PROJECT: Warnings
-------------------------------------------------------------------------
DESCRIPTION:
=========================================================================
%}


clearvars -except projectDir projectDirData fromMainWarningsPaper
clc;close all;fclose('all');feature('DefaultCharacterSet','UTF-8');

if(not(exist('projectDir','var')==1&&exist('projectDirData','var')==1&&exist('fromMainWarningsPaper','var')==1&&fromMainWarningsPaper))
    pcName=char(java.lang.System.getProperty('user.name'));
    if(strcmp(pcName,'felipe'))
        % PC Felipe
        myDir='/Users/felipe/Dropbox/';
        projectDir=[myDir,'git/warnings/'];
        projectDirData=[myDir,'projects/warnings/'];
        addpath(genpath([myDir,'/myMatlabFunctions/']));
    end
end
compileLatexTable=false;


%%
dirPlots=[projectDir,'/paper/figuresCL/RDs/'];
dirTable=[projectDir,'/paper/tablesCL/'];

warning('Alternative dir plots')
dirPlots=[projectDir,'/paper/figuresCL/RDs_aux/'];
dirPlotsR=[projectDir,'/paper/figuresCL/RDs_multBandwidth/'];

dirData=[projectDirData,'/data/chile/'];


tableBandwidthComparison=true; % Make a table comparing estimates with different bandwidths.
posBandwidth=1; % This is for the bandwidht robust excercise. Indicates the subpopulation that we are using the bandwidth table
plotAndSaveRobust=true; % Plots of robust bandwidht RD

plotRDs=true;
titlePlotRD=false;
plotWithZoom=false;
savePlotRDs=true;
bandwidthTable='local'; % 'full','local','localAuto';

switch bandwidthTable
    case 'full'
        bwnote='Full bandwidth was used (-0.28,+0.68) with 2nd order (local) polynomial fit. ';
    case 'local'
        bwnote='Bandwidth was set to +-0.1 with default local polynomial fit (1st order). ';
    case 'localAuto'
        bwnote='Bandwidth was set to automatic (bwselect with default options) with default polynomial fit (1st order). ';
end


% For making subfigures with plots:
maxHorzPlots=3;
maxVertPlots=4;

%% Load data:
anho=0; %0: pooled
if(anho==0)
    anhoStr='';
else
    anhoStr=sprintf('%i',anho);
end


notIn2020={}; % {'assignedToPref','acceptOffer','declineOffer','preferenciaAsign_1','placedInAnyPreferenceIniEnd','placedInAddedIniEnd','placedInAddedNoWaitlistIniEnd','placedInOldIniEnd','enrolledInAssigned','enrolledInAddedIniEnd','enrolledInAddedNoWaitlistIniEnd','enrolledInAddedNoRiskIniEnd'};

if(anho==0)
    
    d1=load([dirData,'2018/inputRD'],'dataRD');
    d2=load([dirData,'2019/inputRD'],'dataRD');
    d3=load([dirData,'2020/inputRD'],'dataRD');
    % Variables that are not available for 2020 yet:
    fillWithNan=notIn2020;
    %fillWithNan=
    for v=1:length(fillWithNan)
        d3.dataRD.(fillWithNan{v})=nan(height(d3.dataRD),1);
    end
    
    
    % Universe that is comparable between 2018-2020
    entryGrades=[-1 0 1 7 9];
    sirven1=ismember(d1.dataRD.grade,entryGrades)&not(d1.dataRD.cod_reg==13);
    sirven2=ismember(d2.dataRD.grade,entryGrades)&not(d2.dataRD.cod_reg==13);
    sirven3=ismember(d3.dataRD.grade,entryGrades)&not(d3.dataRD.cod_reg==13);
    
    d1.dataRD.comparable=sirven1;
    d2.dataRD.comparable=sirven2;
    d3.dataRD.comparable=sirven3;
    
    d1.dataRD.anho=2018*ones(height(d1.dataRD),1);
    d2.dataRD.anho=2019*ones(height(d2.dataRD),1);
    d3.dataRD.anho=2020*ones(height(d3.dataRD),1);
    
    % Keep vars in common
    
    incommon=intersect(intersect(d1.dataRD.Properties.VariableNames,d2.dataRD.Properties.VariableNames),d3.dataRD.Properties.VariableNames);
    incommon=incommon(not(strcmp(incommon,'mrun')));
    dataRD=[d1.dataRD(:,incommon);d2.dataRD(:,incommon);d3.dataRD(:,incommon)];
    
    
    
else
    load([dirData,anhoStr,'/inputRD'],'dataRD')
end

dataRD.rural=dataRD.newMarket<0;

% As asignment at .3 is control, just move (if any) in .3 to epsilon to the
% left:
dataRD.riskPopup(dataRD.riskPopup==.3)=.29999999;

% Avoid mass points at extremes:
dataRD.restrAll=dataRD.pobPopup==1&dataRD.riskPopup>.02&dataRD.riskPopup<.98;

% No need to be more specific, and put "Predicted placement risk of 1st
% attempt"
dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
dataRD.Properties.VariableDescriptions{'treatedPopup'}='Treated pop-up in first attempt';


if(false)
    %% % Popups of total population
    n_apps=[274990,483070,454415];
    
    fprintf('2018 In data: %.1f, in popup population: %.1f\n',height(d1.dataRD)/n_apps(1)*100,sum(d1.dataRD.pobPopup)/n_apps(1)*100);
    fprintf('2019 In data: %.1f, in popup population: %.1f\n',height(d2.dataRD)/n_apps(2)*100,sum(d2.dataRD.pobPopup)/n_apps(2)*100);
    fprintf('2020 In data: %.1f, in popup population: %.1f\n',height(d3.dataRD)/n_apps(3)*100,sum(d3.dataRD.pobPopup)/n_apps(3)*100);
    
    fprintf('\n ALL  In data: %.1f, in popup population: %.1f\n',height(dataRD)/sum(n_apps)*100,sum(dataRD.pobPopup)/sum(n_apps)*100);
    
end





%% RD

% Defines subpopulation in which I want to calculate RDs for outcomes
% defined soon.

if(anho==0) % This calculates estimates pooling all years
    
    
    subpobs=cell(1,4);
    s=1;
    subpobs(s,:)={'Pooled',    dataRD.restrAll,'pooled','addAnyIniEnd'};s=s+1;
    %subpobs(s,:)={'Grade <=8',   dataRD.restrAll&dataRD.grade<9,'basica',''};s=s+1;
    %subpobs(s,:)={'Grade >8',      dataRD.restrAll&dataRD.grade>8,'media',''};s=s+1;
    %subpobs(s,:)={'2020',    dataRD.restrAll&dataRD.anho==2020,'2020',''};s=s+1;
    posAll=1;
    
    
    
    
    
    
else
    error('aca')
end
dataRD.orden=(1:height(dataRD))';

% Definition of outcomes for RDs. 3rd row defines the beginning of a panel
% in the table
% 4th row

if(ismember('DRiskLExpostIniEnd',dataRD.Properties.VariableNames))
    error('Borrar esto');
else
    dataRD.DRiskExpostIniEnd=dataRD.riskExpostEnd-dataRD.riskExpostIni;
end

% New: enrolled conditional on placed (short name because of stata)
dataRD.enrolledInAssignedCOA=double(dataRD.enrolledInAssigned);
dataRD.enrolledInAssignedCOA(dataRD.assignedToPref==0)=nan;
dataRD.placedInAnyPreferenceNWIE=dataRD.placedInAnyPreferenceNoWaitlistIniEnd;
dataRD.placedInAnyPreferenceNoRiskIE=dataRD.placedInAnyPreferenceNoRiskIniEnd;

dataRD.enrolled=not(isnan(dataRD.rbdEnrolled));

%% Extra variables

% Outcomes at individual level:
% data2018=readtable([myDir,'/Mineduc/modelacion/riesgo/data/aggregated/studentsXs2018.csv']);
% data2019=readtable([myDir,'/Mineduc/modelacion/riesgo/data/aggregated/studentsXs2019.csv']);
% data2018.anho=scalarForTable(2018,data2018);
% data2019.anho=scalarForTable(2019,data2019);
% conversion2018=readtable([myDir,'/Mineduc/dataCompilada/Aux/_csv/conversionMRun2018_1.csv'],'delimiter',',');
% conversion2019=readtable([myDir,'/Mineduc/dataCompilada/Aux/_csv/conversionMRun2019_1.csv'],'delimiter',',');
% data2018=outerjoin(data2018,conversion2018,'keys',{'mrun'},'mergeKeys',true,'type','left','rightVariables',{'id_postulante'});
% data2019=outerjoin(data2019,conversion2019,'keys',{'mrun'},'mergeKeys',true,'type','left','rightVariables',{'id_postulante'});
% dataInd=[data2018;data2019];
% 
% assert(mean(ismember(dataInd(:,{'anho','id_postulante'}),dataRD(:,{'anho','id_postulante'}),'rows'))>.94)
% dataRD=outerjoin(dataRD,dataInd,'keys',{'id_postulante','anho'},'mergeKeys',true,'type','left');
% 




dataRD.valueAddedEnrolled(dataRD.grade>8)=nan;
dataRD.hasValueAddedEnrolled=double(not(isnan(dataRD.valueAddedEnrolled)));
dataRD.hasValueAddedEnrolled(dataRD.grade>8)=nan;

dataRD.valueAddedSep=dataRD.valueAddedEnrolled;
dataRD.valueAddedSep(not(dataRD.prioritario))=nan;
dataRD.valueAddedPlaced=dataRD.valueAddedEnrolled;
dataRD.valueAddedPlaced(not(dataRD.placedInAnyPreferenceIniEnd))=nan;
dataRD.valueAddedNotPlaced=dataRD.valueAddedEnrolled;
dataRD.valueAddedNotPlaced((dataRD.placedInAnyPreferenceIniEnd))=nan;


% convCodCursoCodNivel=readtable([myDir,'Mineduc/modelacion/riesgo/data/2020/otherData/relations/eqCodCursoACodNivelFake2019.csv'],'delimiter',',');
% dataLlenos2018=readtable(sprintf('%s/Mineduc/dataCompilada/Aux/_csv/llenos2018.csv',myDir));
% dataLlenos2018.anho=scalarForTable(2018,dataLlenos2018);
% dataLlenos2019=readtable(sprintf('%s/Mineduc/dataCompilada/Aux/_csv/llenos2019.csv',myDir));
% dataLlenos2019.anho=scalarForTable(2019,dataLlenos2019);
% dataLlenos2020=readtable(sprintf('%s/Mineduc/dataCompilada/Aux/_csv/llenos2020.csv',myDir));
% dataLlenos2020.anho=scalarForTable(2020,dataLlenos2020);
% dataLlenos=[dataLlenos2018;dataLlenos2019;dataLlenos2020];
% 
% dataLlenos=outerjoin(dataLlenos,convCodCursoCodNivel,'keys',{'cod_curso'},'mergeKeys',true,'type','left','rightVariables',{'cod_nivel'});
% dataLlenos=stataCollapse({'rbd','cod_nivel','anho'},dataLlenos,'conLista','max');
% dataLlenos=renamevars(dataLlenos,{'rbd','cod_nivel'},{'rbdAsign_1','grade'});
% dataRD=outerjoin(dataRD,dataLlenos,'keys',{'rbdAsign_1','grade','anho'},'mergeKeys',true,'type','left','rightVariables',{'conLista'});
% dataRD=renamevars(dataRD,{'conLista'},{'conListaPlaced'});




dataRD.valueAddedEnrolledCongPlaced=dataRD.valueAddedEnrolled;
dataRD.valueAddedEnrolledCongPlaced(not(dataRD.placedInAnyPreferenceIniEnd))=nan;
dataRD.valueAddedEnrolledCongPlaced(dataRD.conListaPlaced==0)=nan;

dataRD.valueAddedEnrolledUncongPlaced=dataRD.valueAddedEnrolled;
dataRD.valueAddedEnrolledUncongPlaced(not(dataRD.placedInAnyPreferenceIniEnd))=nan;
dataRD.valueAddedEnrolledUncongPlaced(dataRD.conListaPlaced==1)=nan;

% Simce "by grade" Esto viene de "readBaseSimceParaAnalisis.do"
simces=readtable([projectDirData,'/data/chile/auxiliar/simcePerGradeParchado2018.csv'],'delimiter',',');
simces.rbdEnrolled=simces.rbd;
simces.grade=simces.cod_nivel;
simces.simceMat=simces.matParchado/50;
simces.simceLang=simces.langParchado/50;
simces.simcePromedio=(simces.simceMat+simces.simceLang)/2;
dataRD=outerjoin(dataRD,simces,'keys',{'rbdEnrolled','grade'},'mergeKeys',true,'type','left','rightVariables',{'simceMat','simceLang','simcePromedio'});

% 2020 survey data

% From "Mineduc/encuestas/riesgo2020/doFiles/readCleanData.do"
survey=readtable([dirData,'2020/dataEncuestaMail.csv']);
survey.withSurvey=true(height(survey),1);
survey.anho=scalarForTable(2020,survey);
dataRD=outerjoin(dataRD,survey,'keys',{'id_postulante','anho'},'mergeKeys',true,'type','left');

% Check that is not in 100s:
assert(max(dataRD.probPlacedAny)==1)
assert(max(dataRD.probPlaced1st)==1)

dataRD.optimismBiasAny=dataRD.probPlacedAny-(1-dataRD.riskExpostEnd);
dataRD.optimismBias1st=dataRD.probPlaced1st-(1-dataRD.riskExpost_1stPrefEnd);

%%
dataGeorefEEs=readtable(sprintf('%s/data/chile/auxiliar/compilacionInfoGeoref_wNewMarket.csv',projectDirData),'delimiter',',');
dataGeorefEEs=dataGeorefEEs(~isnan(dataGeorefEEs.latitud),:);
dataGeorefEEs.rbdEnrolled=dataGeorefEEs.rbd;
dataGeorefEEs.cod_regEnrolled=dataGeorefEEs.cod_reg;



dataRD=outerjoin(dataRD,dataGeorefEEs,'keys',{'rbdEnrolled'},'type','left','mergekeys',true,'rightVariables',{'cod_regEnrolled','latitud','longitud'});
dataRD.Properties.VariableNames{'latitud'}='latEnrolled';
dataRD.Properties.VariableNames{'longitud'}='lonEnrolled';


dataRD.distanceEnrolled=lldistkm4([dataRD.latEnrolled,dataRD.lonEnrolled],[dataRD.latHome,dataRD.lonHome]);
differentRegion=dataRD.cod_reg~=dataRD.cod_regEnrolled&~isnan(dataRD.cod_regEnrolled);
raros=dataRD.distanceEnrolled>50;

dataRD.distanceEnrolled(differentRegion&raros)=nan;
assert(mean(not(isnan(dataRD.distanceEnrolled)))>.4)



%%
%notIn2020=[notIn2020,'enrolledInAssignedCOA'];




% 5th column defines if IV is calculated or not.
outcomes={'addAnyIniEnd','Add any','a','A. First stage and enrollment',0;...
    'enrolled','Enrolled','e','',0;...
        'hasValueAddedEnrolled','Have value added measure|grade<=8','hva','',0;...
    %'deltaProbPlacedNoRiskIniEnd','$\Delta$ prob. placed to zero risk','dpnr','',1;...
    %'deltaProbPlacedNoWaitlistsIniEnd','$\Delta$ prob. placed to uncongested','dpw','',1;...
    %     'schoolsAddedIniEnd','Schools Added','sa','',1;...
    %     'riskLastDayIni','True Risk initial attempt','ri','';...
    %     'riskLastDayEnd','True Risk final attempt','rf','';...
    %     'DRiskLastDayIniEnd','$\Delta$ True Risk','dr','';...
    %     'riskExpostIni','True Risk initial attempt','riEP','';...
    %     'riskExpostEnd','True Risk final attempt','rfEF','';...
    %     'DRiskExpostIniEnd','$\Delta$ Risk','drEP','',1;...
    %     'numberOfNewUncongestedIniEnd','Uncongested schools Added','saUncong','';...
    %     'addAsFirstIniEnd','Add as first','af','',1;...
    %     'addInBetweenIniEnd','Add to middle','ab','',1;...
    %     'addAsLastIniEnd','Add as last','al','',1;...
    %     'deleteAnyIniEnd','Drop any','da','',1;...
    %     'changeOrigOrderIniEnd','Re-order','ro','',1;...
    %     'riskLastDayEnd','Risk final portfolio','rf','';...
    %     'placedInAnyPreferenceIniEnd','Placed to preference','ap','C. Choice outcome',1;...
    %     'placedInAddedIniEnd','Placed to added preference','apad','';...
    %     'placedInOldIniEnd','Placed to old preference','apold','';...
    %     'enrolledInAssigned','Enrolled in placed','ea','',1;...
    %     'enrolledInAssignedCOA','Enrolled in placed$|$placed','eaadcp','',1;...
    %     'declineOffer','Decline offer','do','',1;...
    %     'enrolledInAddedIniEnd','Enrolled in added','eaad','';...
    %     'addAnyUncongestedIniEnd','Add any uncongested','aaUncong','D. Congestion-related outcomes',1;...
    %     'placedInAnyPreferenceNWIE','Placed to uncong.','apUncong','',1;...
    %     'placedInAddedNoWaitlistIniEnd','Placed to uncong. added pref.','apadUncong','',1;...
    %     'enrolledInAddedNoWaitlistIniEnd','Enrolled in uncong. added','eaadUncong','',1;...
    %     'addAnyNoRiskIniEnd','Add any zero risk','aaNR','D2. Risk-related outcomes',1;...
    %     'placedInAnyPreferenceNoRiskIE','Placed to zero risk','apNR','',1;...
    %     'placedInAddedNoRiskIniEnd','Placed to zero risk added pref.','apadNR','',1;...
    %     'enrolledInAddedNoRiskIniEnd','Enrolled in zero risk added','eaadNR','',1;...
    %     'acceptOffer','Accept offer','ao','';...
    %     'preferenciaAsign_1','Order assigned','oa','';...
    'distanceEnrolled','Distance (km)','dist','B. Attributes of enrolled school',1;...
    %'simceMat','SIMCE Math (score/50)','simce','',0;...
    'valueAddedEnrolled','Value added|grade<=8','va','',1;...
    %'valueAddedPlaced','Value Added$|$placed','vap','',1;...
    %'valueAddedNotPlaced','Value Added$|$not placed','vanp','',1;...
    %'valueAddedEnrolledCongPlaced','Value Added$|$placed to congested','vac','',1;...
    %'valueAddedEnrolledUncongPlaced','Value Added$|$placed to uncongested','vanc','',1;...
    'PerTeacherSpendingEnrolled','Per teacher spending (1000USD)','teachSp','',1;...
    'PerCapSpendingEnrolled','Per student spending (1000USD)','studSp','',1;...
    %'TeacherMath_WeightedAveEnrolled','Teachers'' math score','tms','',1;...
    %'Principal_Score_MathEnrolled','Principal''s score','prs','',1;...
    %'simceMat','Mathematics standarized test','sms','',1;...
    %'simceLang','Language standarized test','sls','',1;...
    %'simcePromedio','Average test score math/reading','aveTest','',1;...
    'hasFeeEnrolled','With copayment fee','hasFee','',1;...
    'feeEnrolled','School monthly fee (USD)','fee','',1;...
    'meanSepEnrolled','Share of vulnerable students','shareVul','',1;...
    %'sizeEnrolled','Total enrollment','se','',1;...
    'sizePerNivelEnrolled','Total enrollment per grade','totalEnrollG','',1;...
    %'privateSubsEnrolled','Subsidized private','pse','',1;...
    %'privateNoSubsEnrolled','Non-subsidized private','pnse','',1;...
    %'conConvenioSepEnrolled','SEP agreement','sepa','',1;...
    %'prom_gral','GPA','gpa','Student level data',1;...
    %'propAsistencia_rend','Attendance 1','at1','',1;...
    %'propAsistencia_asist','Attendance 2','at2','',1;...
    %'valueAddedPlaced','Value Added$|$Placed in preference','vaP','',0;...
    %'valueAddedNotPlaced','Value Added$|$Not placed in preference','vaNP','',0;...
    };

%outcomes=[outcomes;outVA];

% dataRD.highSchool=dataRD.grade>8;
% dataRD.withVA=not(ismissing(dataRD.valueAddedEnrolled));
% 
% % 5th column defines if IV is calculated or not.
% outcomes={'addAnyIniEnd','Add any','a','A. First stage and enrollment',0;...
%     'enrolled','Enrolled','e','',0;...
%     'withVA','With VA measure','wva','',0;...
%     'valueAddedEnrolled','Value Added','va','',1;...
%     };
notIn2020={'prom_gral','propAsistencia_rend','propAsistencia_asist'};

if(anho==2020)
    outcomes=outcomes(not(ismember(outcomes(:,1),notIn2020)),:);
end

dataRD.feeEnrolled=dataRD.feeEnrolled/750;
dataRD.hasFeeEnrolled=double(dataRD.feeEnrolled>0);
dataRD.hasFeeEnrolled(isnan(dataRD.rbdEnrolled))=nan;

dataRD.Properties.VariableDescriptions{'riskPopup'}='Predicted placement risk';
for o=1:size(outcomes,1)
    
    dataRD.Properties.VariableDescriptions{outcomes{o,1}}=outcomes{o,2};
end


assert(allunique(outcomes(:,1)))
assert(allunique(outcomes(:,2)))
assert(allunique(outcomes(:,3)))


% Sort dataRD as it was after merging new vars
dataRD=sortrows(dataRD,'orden');
data=dataRD(:,[{'riskPopup'},outcomes(:,1)']);
inAnyPob=false(height(data),1);




%% Generate stata commands
commands=cell(size(subpobs,1)*size(outcomes,1),2);
counter=1;

for p=1:size(subpobs,1)
    
    data.(sprintf('subpop%i',p))=subpobs{p,2};
    inAnyPob=inAnyPob|subpobs{p,2};
    
    for o=1:size(outcomes,1)
        
        % Avoid computing outcomes that are not available:
        noCalcular=(strcmp(subpobs{p,3},'2020')|strcmp(subpobs{p,3},'2020comp'))&ismember(outcomes{o,1},notIn2020);
        
        %noCalcular=noCalcular|not(strcmp(subpobs{p,3},'2020'))&ismember(outcomes{o,1},surveyOutcomes);
        
        
        
        if(not(noCalcular))
            
            
            % Full bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'full')||plotRDs)
                commands(counter,:)={sprintf('%s_%i_full',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, c(.3) p(2) h(.28 .68)',outcomes{o,1},p)};counter=counter+1;
            end
            % Full loccal fixed bandwidth
            if(tableBandwidthComparison||strcmp(bandwidthTable,'local'))
                commands(counter,:)={sprintf('%s_%i_local',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, c(.3)  h(.1 .1)',outcomes{o,1},p)};counter=counter+1;
            end
            % Full local auto bandwidht
            if(tableBandwidthComparison||strcmp(bandwidthTable,'localAuto'))
                commands(counter,:)={sprintf('%s_%i_localAuto',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, c(.3) ',outcomes{o,1},p)};counter=counter+1;
            end
            
            % If with IV:
            if(outcomes{o,5}==1&&~isempty(subpobs{p,4}))
                
                endogenousVar=subpobs{p,4};
                
                % Full bandwidth
                if(strcmp(bandwidthTable,'full'))
                    commands(counter,:)={sprintf('%s_%i_full_iv',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, fuzzy(%s) c(.3) p(2) h(.28 .68)',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
                % Full loccal fixed bandwidth
                if(strcmp(bandwidthTable,'local'))
                    commands(counter,:)={sprintf('%s_%i_local_iv',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, fuzzy(%s) c(.3)  h(.1 .1)',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
                % Full local auto bandwidht
                if(strcmp(bandwidthTable,'localAuto'))
                    commands(counter,:)={sprintf('%s_%i_localAuto_iv',outcomes{o,3},p),sprintf('rdrobust %s riskPopup if subpop%i==1, fuzzy(%s) c(.3) ',outcomes{o,1},p,endogenousVar)};counter=counter+1;
                end
            end
        end
    end
end
commands=commands(not(cellfun(@isempty,commands(:,1),'UniformOutput',true)),:);


if(false)
    % Run specific rd
    dataRD.auxRest=dataRD.restrAll&dataRD.mandatoryToAdd==0;
    a=stataCommand('rdrobust lengthEnd riskPopup if auxRest==1, c(.3) p(1) h(.1 .1)',dataRD)
    stataExpress('rdrobust  DRiskExpostIniEnd riskPopup if auxRest==1, fuzzy(addAnyIniEnd) c(.3) p(1) h(.1 .1)',data)
    
    yvar=outcomes{1,1};
    rdrobust(yvar,'riskPopup',data.subpop1,.3,' p(2) h(.28 .68)',data);
    figure
    stataExpress('rdplot  enrolled riskPopup if subpop1==1,  c(.3)  p(2) h(.28 .68)',data)
end

% Run Stata:

data=data(inAnyPob,:);
res=stataCommand(commands,data);


%% Make plots and table
close all

cantPobs=(size(subpobs,1));
cantRegs=(size(outcomes,1));

matrix_nl=nan(cantRegs,cantPobs);
matrix_nr=nan(cantRegs,cantPobs);
matrix_b=nan(cantRegs,cantPobs);
matrix_se=nan(cantRegs,cantPobs);

matrix_b_iv=nan(cantRegs,cantPobs);
matrix_se_iv=nan(cantRegs,cantPobs);
matrix_nl_iv=nan(cantRegs,cantPobs);
matrix_nr_iv=nan(cantRegs,cantPobs);

matrix_b_l=nan(cantRegs,cantPobs);
matrix_se_l=nan(cantRegs,cantPobs);

matrix_b_r=nan(cantRegs,cantPobs);
matrix_se_r=nan(cantRegs,cantPobs);

matrix_biasCorr_ci=nan(cantRegs,cantPobs,2);



regNames=fieldnames(res);

cellsubfig=cell(maxHorzPlots,ceil(cantRegs*cantPobs/maxHorzPlots)); % Defined transposed to fill it with scalar indexing.
counter=1;

for p=1:cantPobs
    for r=1:cantRegs
        
        regName_i=sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable);
        if(ismember(regName_i,regNames))
            res_i=res.(regName_i);
            
            % Beta
            matrix_b(r,p)=res_i.tau_cl;
            % Sd Beta
            matrix_se(r,p)=res_i.se_tau_cl;
            
            % Beta left
            matrix_b_l(r,p)=res_i.beta_p_l(1);
            % Sd Beta left
            matrix_se_l(r,p)=sqrt(res_i.V_cl_l(1));
            
            % Beta right
            matrix_b_r(r,p)=res_i.beta_p_r(1);
            % Sd Beta right
            matrix_se_r(r,p)=sqrt(res_i.V_cl_r(1));
            
            % Bias corrected Confidence interval (alpha=5%)
            matrix_biasCorr_ci(r,p,1)=res_i.tau_bc-norminv(.975)*res_i.se_tau_rb;
            matrix_biasCorr_ci(r,p,2)=res_i.tau_bc+norminv(.975)*res_i.se_tau_rb;
            
            % N
            matrix_nl(r,p)=res_i.N_h_l;
            matrix_nr(r,p)=res_i.N_h_r;
            
            % If with IV:
            if(outcomes{r,5}==1&&~isempty(subpobs{p,4}))
                regName_i_iv=sprintf('%s_%i_%s_iv',outcomes{r,3},p,bandwidthTable);
                
                res_i_iv=res.(regName_i_iv);
                % Beta
                matrix_b_iv(r,p)=res_i_iv.tau_cl;
                % Sd Beta
                matrix_se_iv(r,p)=res_i_iv.se_tau_cl;
                
                % N
                matrix_nl_iv(r,p)=res_i_iv.N_h_l;
                matrix_nr_iv(r,p)=res_i_iv.N_h_r;
                
            end
            
            if(plotRDs)
                
                
                res_i=res.(sprintf('%s_%i_%s',outcomes{r,3},p,'full'));
                res_i_table=res.(sprintf('%s_%i_%s',outcomes{r,3},p,bandwidthTable));
                
                plotsub=subpobs{p,2};
                
                figure
                plotRD(res_i,dataRD,'subpop',plotsub,'otherPointEstimate',res_i_table);
                
                if(titlePlotRD)
                    title(sprintf('RD %s - %s',outcomes{r,2},subpobs{p,1}))
                    
                end
                
                if(savePlotRDs)
                    
                    % Guardo pal latexSubfigure
                    filePlot=sprintf(sprintf('%s_P_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3}));
                    cellsubfig{counter}=  easyExport([dirPlots,filePlot],'title',filePlot);
                    counter=counter+1;
                end
                
                
                
                if(plotWithZoom)
                    % Plot with zoom
                    figure
                    plotRD(res_i,dataRD,'subpop',plotsub,'newxlim',[.1,.5],'nb',100,'otherPointEstimate',res_i_table);
                    easyExport([dirPlots,sprintf(sprintf('%s_P_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                end
                
                
                
            else
                fprintf('Ojo: reg %s no existe para subpob %s\n',regName_i,subpobs{p,1});
            end
        end
    end
    
    %vector_Nl(1,p)=sum(dataRD.riskPopup<.3&subpobs{p,2});
    %vector_Nr(1,p)=sum(dataRD.riskPopup>.3&subpobs{p,2});
    
end





%% RD main table
matTable=nan(size(matrix_b).*[1 2]+[2 0]);
matTable_sd=nan(size(matrix_b).*[1 2]+[2 0]);

matTable(:,1:2:end)=[matrix_b;...
    matrix_nl(end,:)...
    ;matrix_nr(end,:)];

matTable(:,2:2:end)=[matrix_b_iv;...
    matrix_nl_iv(find(cell2mat(outcomes(:,5)),1,'last'),:)...
    ;matrix_nr_iv(find(cell2mat(outcomes(:,5)),1,'last'),:)];


matTable_sd(:,1:2:end)=[matrix_se;nan(2,size(matrix_se,2))];
matTable_sd(:,2:2:end)=[matrix_se_iv;nan(2,size(matrix_se_iv,2))];

withInfo=not(all(isnan(matTable),1));


header=repmat({''},2,size(matTable,2));
header(1,1:2:end)=subpobs(:,1)';
header(1,2:2:end)=subpobs(:,1)';
header(2,2:2:end)={'IV'};

matTable=matTable(:,withInfo);
matTable_sd=matTable_sd(:,withInfo);
header=header(:,withInfo);
header=header(any(not(cellfun(@isempty,header)),2),:);

% add mean
matTable=[matTable,[matrix_b_l(:,posAll);nan;nan]];
matTable_sd=[matTable_sd,nan(size(matTable_sd,1),1)];
header=[header,{'';'$E[Y|X=0.3^{-}]$'}];

opt=struct;
opt.header=header;
opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=[outcomes(:,2);{'NL';'NR'}];
opt.addColumnNumber=true;
if(anho==0)
    opt.columnafantasma=4;
    
else
    %opt.columnafantasma=[4];
    opt.title=sprintf('%s Sample - Pop-Up effect',anhoStr);
end
opt.filaFantasma=size(matTable,1)-2;
opt.columnaFantasma=1;
paneles=find(not(ismissing(outcomes(:,end-1))));
opt.panel=[num2cell(paneles-1),outcomes(paneles,end-1)];
%opt.alignmentFirstCol={'L{3cm}'};
opt.adjust=true;

opt.positionParameter='H';
opt.file=sprintf('%s%s_tablePopup_enrolledOtucomes',dirTable,anhoStr);

opt.title='RD Estimates of Platform Pop-Up Effects on Enrolled School Outcomes';
opt.label='tab:RDpopup-enrolledOutcomes';

opt.note='	Local linear RD estimates of platform pop-up effects. Computed using triangular kernel with bandwidth 0.1. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}. IV estimates in the second column report instrumental variable specifications where the endogenous regressor is the ``add any school'''' indicator. The third column reports below-cutoff means of the variable listed in the row in the analysis sample. Sample for value added outcomes is restricted to grades eight and below. Reported sample sizes are counts of enrolling students. See section \ref{sec:schoolXRD} for discussion and Online Appendix \ref{sec:schoolattappendix} for detailed variable definitions';
tabla=cell2latex(mat2cellstr(matTable,'precisionDecimal','%.3f','revisarFilas',true),'opts',opt);
if(compileLatexTable)
compileLatex(tabla)
end

%%

if(plotRDs&&savePlotRDs)
    close all
    numOfPages=ceil(size(cellsubfig,2)/maxVertPlots);
    pages=cell(numOfPages,1);
    for p=1:numOfPages
        
        newEnd=min((p*maxVertPlots),size(cellsubfig,2));
        note='';
        %note='Pop-up warning RD effect fits and point estimates by bandwidth for outcomes listed in panel titles. ``Full'''': global quadratic. ``+/- 0.1'''': local linear within 0.1 bandwidth. ``rdbwselect'''': optimal bandwidth selection using \citet{calonico2014robust}.';
        if p==1
            label='fig:appBW';
        else
            label=sprintf('fig:appBW_%i',p);
        end
        label='';
        
        %pages{p}=subfiguresLatex(cellsubfig(:,((p-1)*maxVertPlots+1):newEnd)','erp','','erpf', '','caption','RD plots',...
        %    'scale',.45,'docwidth',1.2,'note',note,'label',label);
        
    end
    fileSubfigure='subfigRDs_other';
    if(compileLatexTable)
    %compileLatex(horzcat(tabla,newline,'\clearpage',newline,pages{:}),'dirTex',dirPlots,'texFile',fileSubfigure);
    end
end



%% Make table of bandwidth comparison

if(tableBandwidthComparison)
    cantRobuts=3;
    labelRobust={'Full','+-0.1','rdbwselect'};
    cantRegs=(size(outcomes,1));
    
    matrix_nl=nan(cantRegs,cantRobuts);
    matrix_nr=nan(cantRegs,cantRobuts);
    matrix_b=nan(cantRegs,cantRobuts);
    matrix_se=nan(cantRegs,cantRobuts);
    
    matrix_bl=nan(cantRegs,cantRobuts);
    matrix_br=nan(cantRegs,cantRobuts);
    
    regNames=fieldnames(res);
    
    
    cellsubfig=cell(maxHorzPlots,ceil(cantRegs/maxHorzPlots)); % Defined transposed to fill it with scalar indexing.
    counter=1;
    for p=posBandwidth
        for r=1:cantRegs
            
            regName_i=sprintf('%s_%i_full',outcomes{r,3},p);
            regName_i_local=sprintf('%s_%i_local',outcomes{r,3},p);
            regName_i_localAuto=sprintf('%s_%i_localAuto',outcomes{r,3},p);
            
            if(ismember(regName_i_local,regNames))
                res_i=cell(cantRobuts,1);
                
                res_i{1}=res.(regName_i);
                res_i{2}=res.(regName_i_local);
                res_i{3}=res.(regName_i_localAuto); % AUTO tiene que ser la ultima!
                
                
                for i=1:cantRobuts
                    % Beta
                    matrix_b(r,i)=res_i{i}.tau_cl;
                    % Sd Beta
                    matrix_se(r,i)=res_i{i}.se_tau_cl;
                    
                    % N
                    matrix_nl(r,i)=res_i{i}.N_h_l;
                    matrix_nr(r,i)=res_i{i}.N_h_r;
                    
                    % N
                    matrix_bl(r,i)=res_i{i}.h_l;
                    matrix_br(r,i)=res_i{i}.h_r;
                    
                end
                
                if(plotAndSaveRobust)
                    
                    plotsub=subpobs{p,2};
                    
                    
                    
                    % Plot alternatives bandwidths:
                                        linestylesAB={'-','--','-.'};
                    colors=linspecer(3);
                    c1=colors(1,:);
                    c2=colors(2,:);
                    c3=colors(3,:);
                    
                    figure
                    % Plot full bandwidth:
                    
                    a1=plotRD(res_i{1},dataRD,'subpop',plotsub,'color',c1,'linestyle',linestylesAB{1} ); hold on;
                    if(a1.posBeta>.5);shift=-.1;else;shift=.1;end
                    a2=plotRD(res_i{2},dataRD,'subpop',plotsub,'withCutoffVerticalLine',0,'withBinsreg',0,'posbeta',a1.posBeta+shift,'color',c2,'linestyle',linestylesAB{2} );  hold on;
                    a3=plotRD(res_i{3},dataRD,'subpop',plotsub,'withCutoffVerticalLine',0,'withBinsreg',0,'posbeta',a1.posBeta+2*shift,'color',c3,'linestyle',linestylesAB{3} );
                    
                    legend([a1.functionLine,a2.functionLine,a3.functionLine],labelRobust,'interpreter','latex')
                    
                    hold off
                    
                    % Guardo pal latexSubfigure
                    cellsubfig{counter}=easyExport([dirPlotsR,sprintf(sprintf('%s_P_%s_%s_bandComparison%s.png',anhoStr,outcomes{r,3},subpobs{p,3}))],'caption',outcomes{r,2});
                    counter=counter+1;
                    
                    
                    close
                    
                    %                 % Plot with zoom
                    %                 figure
                    %                 a=plotRD(res_i,dataRD,plotsub,'newxlim',[.1,.5]);
                    %                 easyExport([dirPlots,sprintf(sprintf('%s_P_%s_%s_withZoom.png',anhoStr,outcomes{r,3},subpobs{p,3}))]);
                    %
                    %                 close
                    
                end
            end
        end
        
        
    end
    
   %%
    if(plotAndSaveRobust)
    for p=1:ceil(size(cellsubfig,2)/maxVertPlots)
        
        newEnd=min((p*maxVertPlots),size(cellsubfig,2));
        note='Pop-up warning RD effect fits and point estimates by bandwidth for outcomes listed in panel titles. ``Full'''': global quadratic. ``+/- 0.1'''': local linear within 0.1 bandwidth. ``rdbwselect'''': optimal bandwidth selection using \citet{calonico2014robust}. See  section \ref{sec:RD} for details.';
        if p==1
            label='fig:appBW-enrolledOutcomes';
        else
            label=sprintf('fig:appBW-enrolledOutcomes-%i',p);
        end
        %subfiguresLatex(cellsubfig(:,((p-1)*maxVertPlots+1):newEnd)','file',[dirPlotsR,'containers/subfigRobustRDs-enrolledOtucomes_',num2str(p)],'erp','figuresCL/RDs_multBandwidth/','erpf', 'figuresCL/RDs_multBandwidth/','caption','Multiple bandwidths RD plots of platform-based pop-up warning effects (outcomes in Table \ref{tab:RDpopup-enrolledOutcomes})',...
        %    'scale',.52,'docwidth',1.2,'note',note,'label',label);
    end
    end
    %%
    
    cellTable=[[mat2cellstr(matrix_b,'precisionDecimal','%.3f','revisarCols',true),mat2cellstr([matrix_bl(:,end),matrix_br(:,end),matrix_nl(:,end),matrix_nr(:,end)],'precisionDecimal','%.2f','revisarCols',true)];...
        mat2cellstr([[matrix_nl(end,1:2);matrix_nr(end,1:2)],nan(2,5)])];
    
    nans=nan(size(cellTable));
    nans(1:size(matrix_se,1),1:size(matrix_se,2))=matrix_se;
    matTable_sd=nans;
    
    
    opt=struct;
    opt.header=[[{'Bandwidth'},labelRobust,{'rdbwselect','rdbwselect','rdbwselect','rdbwselect'}];...
        {'','Estimate','Estimate','Estimate','BW left','BW right','N left','N right'}];
    opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
    opt.primeracolumna=[outcomes(:,2);{'N left';'N right'}];
    %opt.stars=getStars(matTable,matTable_sd);
    opt.addColumnNumber=true;
    opt.columnafantasma=[1 2];
    opt.title='Platform Pop-Up RD Estimates of Enrolled School Outcomes (Table \ref{tab:RDpopup-enrolledOutcomes}) with Alternate Bandwidths';
    opt.filaFantasma=size(outcomes,1);
    opt.positionParameter='H';
    paneles=find(not(ismissing(outcomes(:,4))));
    opt.panel=[num2cell(paneles-1),outcomes(paneles,4)];
    %opt.alignmentFirstCol={'L{3cm}'};
    opt.adjust=true;
    opt.label='tab:popupBW-enrolledOutcomes';
    opt.file=sprintf('%s%s_tablePopupBandwidth-enrolledOutcomes',dirTable,anhoStr);
    opt.note="Local linear and full sample quadratic polynomial RD estimates of pop-up effects from warning pop-up on application platform. Computed using triangular kernel with different bandwidths. ``Full'' bandwidth uses 2nd order polynomial fit, while ``+-0.1'' and ``rdbwselect'' use 1st order (local) polynomial. Heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors reported in parentheses; computed as in \citet{calonico2014robust}.";
    tablaComp=cell2latex(cellTable,'opts',opt);
    if(compileLatexTable)
    compileLatex(tablaComp)
    end
    
end
